%load data
close all;clc;
load '../data/lsLNH.mat';

ls1M = val(:,2);
lls1M  = ls1M;
% lls1M = log(val(:,2));
% rls1M = log(ls1M(2:end)./ls1M(1:end-1));
% [transR, lambda] = boxcox(ls1M);

%%
rlls1M = lls1M(2:end)-lls1M(1:end-1);

%% check ARCH
close all;
figure(4)
plot(rlls1M);

N = size(rlls1M,1);

figure
subplot(2,1,1)
autocorr(rlls1M)
subplot(2,1,2)
parcorr(rlls1M)

[h,p] = lbqtest(rls1M,'Lags',5)

%% Step 2. Check the series for autocorrelation.
e = rlls1M - mean(rlls1M);
figure
subplot(2,1,1)
autocorr(e.^2)
subplot(2,1,2)
parcorr(e.^2)

%% Step 3. Test the significance of the autocorrelations
[h,p] = lbqtest(e.^2,'Lags',[1:10]);
[h' p']
% the value is reject the noncorrelation -> there is GARCH


%% Step 4. Check the series for conditional heteroscedasticity.
figure
subplot(2,1,1)
autocorr(rlls1M.^2)
subplot(2,1,2)
parcorr(rlls1M.^2)


%% Step 5. Test for significant ARCH effects

% [h,p] = archtest(rls1M-mean(rls1M),'lags',1:5)

numLags = 6;
% logL = zeros(numLags,1); % Preallocate fit statistics
res = zeros(numLags+1,numLags);
for h = 0: numLags
for k = 1:numLags
    Mdl = garch(h,k); % Specify garch model
    [~,~,logL] = estimate(Mdl,rlls1M,'display','off'); % Obtain loglikelihood
    [~,res(h+1,k)] = aicbic(logL,h+k,N);
end
end
% end

res

% fitStats = aicbic(logL,[1:numLags]) % Get AIC
% lags = find(min(fitStats)) % Obtain suitable number of lags

%%
% Conduct the ARCH Test at signifiance 5%
[h,pValue,stat,cValue] = archtest(rls1M-mean(rls1M),'lags',lags,'alpha',0.05);
[h,pValue,stat,cValue]

%% Step 6. Specify a conditional mean and variance model.
model = arima('D',1,'ARLags',[],'MALags',[],'Variance',garch(1,1));


% Step 7: Estimate the model parameters without using presample data.
fit = estimate(model,lls1M)



%%

[res,V,LogL] = infer(fit,lls1M);

figure
subplot(2,1,1)
plot(V)
xlim([0,N])
title('Conditional Variance')

subplot(2,1,2)
plot(res./sqrt(V))
xlim([0,N])
title('Standardized Residuals')


%% %% Step 8. Fit a model with a t innovation distribution.
% 
% model.Distribution = 't';
% fitT = estimate(model,rls1M,'Variance0',{'Constant0',0.001});
% 
% %% loop to get the best model using AIC

%%
Npredict = 20;
[Y,YMSE,V1] = forecast(fit,Npredict,'Y0',lls1M,'E0',res,'V0',V);
upper = Y + 1.96*sqrt(YMSE);
lower = Y - 1.96*sqrt(YMSE);
Nstart = 200;

figure
subplot(2,1,1)
plot(lls1M(end-Nstart+1:end),'Color',[.75,.75,.75])
hold on
plot(Nstart+1:Nstart+Npredict,Y,'r','LineWidth',2)
plot(Nstart+1:Nstart+Npredict,[upper,lower],'k--','LineWidth',1.5)
xlim([0,Nstart+Npredict])
title('Forecasted Returns')
hold off
subplot(2,1,2)
plot(V(end-Nstart+1:end),'Color',[.75,.75,.75])
hold on
plot(Nstart+1:Nstart+Npredict,V1,'r','LineWidth',2);
xlim([0,Nstart+Npredict])
title('Forecasted Conditional Variances')
hold off


%% simulate
% rng('default')
% 
% 
% [Y1,E1,V1] = simulate(fit,Npredict,'NumPaths',10,...
%                    'Y0',rls1M,'E0',res,'V0',V);
% 
% Y0 = ls1M(end);
% 
% r = zeros(Npredict+1,10);
% for i = 1: Npredict    
%     r(i+1,:) = exp(log(Y0) + Y1(i,:));    
%     Y0 = r(i+1);
% end     
% r(1,:) = ls1M(end);
%                
% figure
% plot(ls1M,'Color',[.75,.75,.75])
% hold on
% plot(N:N+Npredict,r(:,9))
% xlim([0,N+Npredict])
% title('Simulated Returns')
% hold off

